function sol=estabilidad()
h = .1;
sol = [];
RTc = .05:.05:.85;
RTv = .1:.05:.9;
for [amplificacion, name] = struct(
		'euler',	@(Tc, Tv, h) (1 + h*(Tc - Tv)),
		'rk2',		@(Tc, Tv, h) (1 + h*(Tc - Tv)*(2 + Tc - Tv)/2),
		'rk4',		@(Tc, Tv, h) (1 + (Tc - Tv)*(((h*(Tc-Tv))^3)/24
													+ ((h*(Tc-Tv))^2)/6
													+ (h*(Tc-Tv))/2
													+ 1)))
		s = [];
		m = [];
		for Tc = RTc
				r = [];
				for Tv = RTv
						a = amplificacion(Tc, Tv, h);
						if (Tc != Tv)
								s = [s; Tc, Tv, a];
						end;
						r = [r,a];
				end;
				m = [m; r];
		end;
		sol.(name).FAs = sortrows(s, 3);
		bestconds = sol.(name).FAs(1, :);
		sol.(name).bestconds = struct('Tc', bestconds(1), 'Tv', bestconds(2));
		sol.(name).mesh = m;
		figure();
		meshc(RTv, RTc, m);
		xlabel('Tv');
		ylabel('Tc');
		zlabel('FA');
		print(sprintf('estabilidad_%s.png', name));
end;
end;
